This system w single functional group, “plants.” Carbon biomass as currency. Carbon biomass density as \(P(t)\) gC/m2. Photosynthesis produces new biomass at \(\phi\) gC/m2/day, and plant of mass \(w\) loses mass (mortality and respiration) at \(\delta_p w\) gC/day.
\[\frac{dP}{dt} = \phi - \delta_p P\] So, steady state stock of plants at \[0 = \phi - \delta_p P \Longrightarrow P^* = \phi / \delta_p\]
Add a second group, “herbivores,” to the ecosystem. Carbon biomass density of herbivores is \(H(t)\) gC/m2; respiration & mortality loss at \(\delta_h\) day-1, and linear functional response attack rate \(\alpha_h\) m^2/day/gC.
\[\frac{dP}{dt} = \phi - \delta_p P - \alpha_h PH\] and \[\frac{dH}{dt} = \alpha_h PH - \delta_h H\] Can find two steady states, one with \(H=0\) and \(P\) given by one-level steady state above. The other is: \[P^* = \delta_h / \alpha_h, \quad H^* = \frac{1}{\delta_h}(\phi - \delta_p P^*)\]
This will calculate \(\frac{dH}{dt}\) and \(\frac{dP}{dt}\).
### params for simulations
a_h <- 0.01 ### herbivore attack rate
phi <- 10 ### photosynthesis rate gC/m^2/day
d_h <- 0.5 ### carbon mass losses for herbivores
d_p <- 0.1 ### carbon mass losses for plants
pops_0 <- c(P = 25, H = 5) ### initial carbon densities of herbivore and plant
prms <- c(a_h = a_h, phi = phi, d_h = d_h, d_p = d_p)
t <- seq(0, 100, by = 1)
foodchain2_out <- ode(pops_0, t, foodchain_2level, prms)
foodchain2_df <- foodchain2_out %>%
as.data.frame() %>%
setNames(c('t', 'P', 'H'))### Equilibrium
P_star2 <- d_h / a_h
H_star2 <- 1/d_h * (phi - d_p * P_star2)
### Isoclines
P_iso2 <- function(p) {
h = (phi - d_p * p) / (a_h * p)
}
H_iso2 <- P_star2
### plot
ggplot(data = foodchain2_df, aes(x = P, y = H)) +
ggtheme_plot() +
stat_function(aes(x = P), fun = P_iso2,
linetype = 'dashed', color = 'darkred') +
geom_vline(xintercept = H_iso2,
linetype = 'dashed', color = 'darkblue') +
geom_segment(aes(xend = lead(P), yend = lead(H)),
arrow = arrow(length = unit(0.07, 'inches')),
size = .25) +
geom_point(x = P_star2, y = H_star2, size = 3) +
geom_hline(yintercept = 0) + geom_vline(xintercept = 0) +
ylim(c(0, max(foodchain2_df$H)))This is not in Ecological Dynamics, but based on Primer.
From the equil formulas above, \(P^* = \frac{\delta_h}{\alpha_h}\) = 50, \(H^* = \frac{1}{\delta_h}(\phi - \delta_p P^*)\) = 10
To know how each population growth rate changes in response to changes in the abundance each of the other population, we take the partial derivatives.
\[\mathbf J = \begin{pmatrix} \frac{\partial \dot P}{\partial P} & \frac{\partial \dot P}{\partial H} \\ \frac{\partial \dot H}{\partial P} & \frac{\partial \dot H}{\partial H} \end{pmatrix} = \begin{pmatrix} -\delta_p - \alpha_h H & -\alpha_h P \\ \alpha_h H & \alpha_h P - \delta_h \end{pmatrix}\]
We can replace the \(P\) and \(H\) in the Jacobian with the equilibria found above \(P^*, H^*\).
\[\begin{pmatrix} -\delta_p - \alpha_h \frac{1}{\delta_h}(\phi - \delta_p \frac{\delta_h}{\alpha_h}) & -\alpha_h \frac{\delta_h}{\alpha_h} \\ \alpha_h \frac{1}{\delta_h}(\phi - \delta_p \frac{\delta_h}{\alpha_h}) & \alpha_h \frac{\delta_h}{\alpha_h} - \delta_h \end{pmatrix} = \begin{pmatrix} \frac{-\alpha_h \phi}{\delta_h} & -\delta_h\\ \frac{\alpha_h \phi}{\delta_h} - \delta_p & 0 \end{pmatrix}\]
Recall the Routh-Hurwitz criterion:
\[ J_{11} + J_{22} < 0\] \(J_{11}\) is negative, so the plant population will self regulate, while \(J_{22}\) is zero, so the herbivore population does not have a density-dependence effect on itself. ???
The other part of the Routh-Hurwitz criterion is the condition:
\[ J_{11}J_{22} - J_{12}J_{21} > 0\] In the plant-herbivore context, this suggests that the plant declines due to the herbivore (\(J_{12} < 0\)) and the herbivore increases due to the plant (\(J_{21} > 0\) as long as \(\phi > \delta_p \frac{\delta_h}{\alpha_h}\)). The signs of these elements make their product negative, and help make the above condition true.
Note: the \(\phi > \frac{\delta_p \delta_h}{\alpha_h}\) term indicates: primary production \(\phi\) is high enough to support respiration/mortality losses from the plant population at the density needed to bring the herbivore production and loss into balance, or alternately, there’s only a finite herbivore standing stock if the plant standing stock in the absennce of herbivory is larger than that required to sustain steady-state herbivore standing stock.
# We can perform eigenanalysis given the parameters above.
Jac <- matrix(c(-a_h * phi / d_h, -d_h, a_h * phi / d_h - d_p, 0),
byrow = TRUE, nr = 2)
eig <- round(eigen(Jac)[["values"]],3)If we performed eigenanalysis on the above Jacobian matrix, we would find that the eigenvalues are complex (-0.1+0.2i, -0.1-0.2i), meaning that the populations will oscillate or cycle, with period \(\frac{2\pi}{\omega}\). * If the real parts are zero, the system will exhibit neutral stability * If the real parts are negative, the system will converge toward the equilibrium. * If the real parts are positive, the system will approach a stable oscillation away from the equilibrium.
Add a third group, “consumers,” which only eat herbivores, to the ecosystem. Carbon biomass density of consumers is \(C(t)\) gC/m2; respiration & mortality loss at \(\delta_c\) day-1, and linear functional response attack rate \(\alpha_c\) m^2/day/gC. The system of equations:
\[\frac{dP}{dt} = \phi - \delta_p P - \alpha_h PH\] \[\frac{dH}{dt} = \alpha_h PH - \delta_h H - \alpha_c H C\] \[\frac{dC}{dt} = \alpha_c HC - \delta_c C\] Can find three steady states:
This will calculate \(\frac{dC}{dt}, \frac{dH}{dt}, \frac{dP}{dt}\).
To have sensible equilibria, the following needs to hold true: \[\phi \geq \delta_p \left(\frac{\delta_h}{\alpha_h}\right) + \delta_h \left(\frac{\delta_c}{\alpha_c} \right)\]
### params for simulations: what are some sensible numbers?
a_h <- 0.005 ### herbivore attack rate
a_c <- 0.015 ### consumer attack rate
phi <- 150 ### photosynthesis rate gC/m^2/day
d_h <- 0.90 ### carbon mass losses for herbivores
d_p <- 0.20 ### carbon mass losses for plants
d_c <- 0.93 ### carbon mass losses for consumers
phi_floor <- d_p * d_h / a_h + d_h * d_c / a_c
if(phi < phi_floor) {
stop('phi value is too low')
}
### Equilibrium
H_star3 <- d_c / a_c
P_star3 <- phi / (d_p + a_h * H_star3)
C_star3 <- (phi - d_p * P_star3 - d_h * H_star3) / d_c
if(any(c(P_star3 < 0, H_star3 < 0, C_star3 < 0))) {
stop('Param values lead to non-sensible equilibria')
}pops_0 <- c(P = 100, H = 100, C = 5)
### initial carbon densities of herbivore, plant, and consumer
prms <- c(a_h = a_h, a_c = a_c, phi = phi,
d_h = d_h, d_p = d_p, d_c = d_c)
t <- seq(0, 100, by = 0.1)
foodchain3_out <- ode(pops_0, t, foodchain_3level, prms)
foodchain3_df <- foodchain3_out %>%
as.data.frame() %>%
setNames(c('t', 'P', 'H', 'C'))equil_df <- data.frame(P = P_star3, H = H_star3, C = C_star3) %>%
gather(level, pop) %>%
mutate(level = fct_inorder(level))
pop_v_time_df <- foodchain3_df %>%
gather(level, pop, -t) %>%
mutate(level = fct_inorder(level))
ggplot(pop_v_time_df, aes(x = t, y = pop, color = level)) +
ggtheme_plot() +
geom_line() +
geom_hline(data = equil_df, aes(yintercept = pop, color = level),
linetype = 'dashed', size = .25) +
geom_hline(yintercept = 0) +
scale_color_manual(values = c('darkgreen', 'blue', 'darkred'))From the equil formulas above:
To know how each population growth rate changes in response to changes in the abundance each of the other population, we take the partial derivatives.
\[\mathbf J = \begin{pmatrix} \frac{\partial \dot P}{\partial P} & \frac{\partial \dot P}{\partial H} & \frac{\partial \dot P}{\partial C} \\ \frac{\partial \dot H}{\partial P} & \frac{\partial \dot H}{\partial H} & \frac{\partial \dot H}{\partial C} \\ \frac{\partial \dot C}{\partial P} & \frac{\partial \dot C}{\partial H} & \frac{\partial \dot C}{\partial C} \end{pmatrix} = \begin{pmatrix} -\delta_p - \alpha_h H & -\alpha_h P & 0\\ \alpha_h H & \alpha_h P - \delta_h - \alpha_c C & -\alpha_c H\\ 0 & \alpha_c C & \alpha_c H - \delta_c \end{pmatrix}\]
We can replace the \(P, H, C\) in the Jacobian with the equilibria found above \(P^*, H^*, C^*\).
\[\mathbf J = \begin{pmatrix} -\delta_p - \frac{\alpha_h \delta_c}{\alpha_c} & - \frac{\alpha_h \phi}{\delta_p + \frac{\alpha_h \delta_c}{\alpha_c}} & 0\\ \frac{\alpha_h \delta_c}{\alpha_c} & \frac{\alpha_h \phi}{\delta_p + \alpha_h \frac{\delta_c}{\alpha_c}} - \delta_h - \frac{\alpha_c }{\delta_c}(\phi - \frac{\delta_p \phi}{\delta_p + \frac{\alpha_h \delta_c}{\alpha_c}} - \frac{\delta_h\delta_c}{\alpha_c}) & -\frac{\alpha_c \delta_c}{\alpha_c}\\ 0 & \frac{\alpha_c}{\delta_c}(\phi - \frac{\delta_p \phi}{\delta_p + \frac{\alpha_h \delta_c}{\alpha_c}} - \frac{\delta_h \delta_c}{\alpha_c}) & \frac{\alpha_c \delta_c}{\alpha_c} - \delta_c \end{pmatrix}\] \[\Longrightarrow \mathbf J = \begin{pmatrix} -\delta_p - \frac{\alpha_h\delta_c}{\alpha_c} & - \frac{\alpha_h \phi}{\delta_p + \frac{\alpha_h \delta_c}{\alpha_c}} & 0\\ \frac{\alpha_h \delta_c}{\alpha_c} & 0 & -\delta_c\\ 0 & \frac{\alpha_h \phi}{\delta_p + \frac{\alpha_h \delta_c}{\alpha_c}} - \delta_h & 0 \end{pmatrix}\]
Recall the Routh-Hurwitz criterion, assuming we can update simply for three species?:
\[ J_{11} + J_{22} + J_{33}< 0\] \(J_{11}\) is negative, so the plant population will self regulate, while \(J_{22}\) and \(J_{33}\) are zero, so the herbivore and consumer populations do not have a density-dependence effect on themselves. ???
The other part of the Routh-Hurwitz criterion is the condition (how to update for three species?):
\[ J_{11}J_{22} - J_{12}J_{21} > 0\]
# We can perform eigenanalysis given the parameters above.
Jac3 <- matrix(
c(-d_p - a_h*d_c / a_c, - (a_h*phi) * (d_p + a_h*d_c / a_c), 0,
a_h*d_c / a_c, 0, -d_c,
0, (a_h * phi) * (d_p + (a_h*d_c) / a_c) - d_h, 0),
byrow = TRUE, nr = 3)
eig3 <- round(eigen(Jac3)[["values"]],3)When we performed eigenanalysis on the above Jacobian matrix, we find that the eigenvalues are complex (0.644+0i, -0.577+0.219i, -0.577-0.219i), meaning that the populations will oscillate or cycle. What does it mean that one eigenvalue is not complex? * If the real parts are zero, the system will exhibit neutral stability * If the real parts are negative, the system will converge toward the equilibrium. * If the real parts are positive, the system will approach a stable oscillation away from the equilibrium.
Add a fourth group, “top predators,” which only eat consumers, to the ecosystem. Carbon biomass density of top predators is \(T(t)\) gC/m2; respiration & mortality loss at \(\delta_t\) day-1, and linear functional response attack rate \(\alpha_t\) m^2/day/gC. The system of equations:
\[\frac{dP}{dt} = \phi - \delta_p P - \alpha_h PH\] \[\frac{dH}{dt} = \alpha_h PH - \delta_h H - \alpha_c H C\] \[\frac{dC}{dt} = \alpha_c HC - \delta_c C - \alpha_c CT\] \[\frac{dT}{dt} = \alpha_t CT - \delta_t T\] Can find four steady states:
This will calculate \(\frac{dC}{dt}, \frac{dH}{dt}, \frac{dP}{dt}\).
foodchain_4level <- function(t, y, params) {
P <- y[1]
H <- y[2]
C <- y[3]
TT <- y[4]
with(as.list(params), {
dP_dt <- phi - d_p * P - a_h * P * H
dH_dt <- a_h * P * H - d_h * H - a_c * H * C
dC_dt <- a_c * H * C - d_c * C - a_t * C * TT
dT_dt <- a_t * C * TT - d_t * TT
return(list(c(dP_dt, dH_dt, dC_dt, dT_dt)))
})
}To have sensible equilibria, the following needs to hold true: ??? update \[\phi \geq \delta_p \left(\frac{\delta_h}{\alpha_h}\right) + \delta_h \left(\frac{\delta_c}{\alpha_c} \right)\]
### params for simulations: what are some sensible numbers?
a_h <- 0.015 ### herbivore attack rate
a_c <- 0.025 ### consumer attack rate
a_t <- 0.035 ### top predator attack rate
phi <- 150 ### photosynthesis rate gC/m^2/day
d_h <- 0.90 ### carbon mass losses for herbivores
d_p <- 0.25 ### carbon mass losses for plants
d_c <- 0.90 ### carbon mass losses for consumers
d_t <- 0.90 ### carbon mass losses for top predators
phi_floor <- d_p * d_h / a_h + d_h * d_c / a_c
if(phi < phi_floor) {
stop('phi value is too low')
}
### Equilibrium
C_star4 <- d_t / a_t
P_star4 <- (d_h + a_c * C_star4) / a_h
H_star4 <- (phi - d_p * P_star4) / (a_h * P_star4)
T_star4 <- (a_c * H_star4 - d_c) / (a_t)
if(any(c(P_star4 < 0, H_star4 < 0, C_star4 < 0, T_star4 < 0))) {
stop('Param values lead to non-sensible equilibria')
}pops_0 <- c(P = 100, H = 100, C = 20, TT = 10)
### initial carbon densities of herbivore, plant, and consumer
prms4 <- c(a_h = a_h, a_c = a_c, a_t = a_t, phi = phi,
d_h = d_h, d_p = d_p, d_c = d_c, d_t = d_t)
t <- seq(0, 100, by = 0.1)
foodchain4_out <- ode(pops_0, t, foodchain_4level, prms4)
foodchain4_df <- foodchain4_out %>%
as.data.frame() %>%
setNames(c('t', 'P', 'H', 'C', 'T'))equil_df <- data.frame(P = P_star4, H = H_star4, C = C_star4, T = T_star4) %>%
gather(level, pop) %>%
mutate(level = fct_inorder(level))
pop_v_time_df <- foodchain4_df %>%
gather(level, pop, -t) %>%
mutate(level = fct_inorder(level))
ggplot(pop_v_time_df, aes(x = t, y = pop, color = level)) +
ggtheme_plot() +
geom_line() +
geom_hline(data = equil_df, aes(yintercept = pop, color = level),
linetype = 'dashed', size = .25) +
geom_hline(yintercept = 0) +
scale_color_manual(values = c('darkgreen', 'blue', 'red', 'darkred'))Modify plant carbon biomass to grow logistically to carrying capacity \(K\):
\[\frac{dP}{dt} = rP \left( 1 - \frac{P}{K} \right)\]
Add a second group, “herbivores,” to the ecosystem, following assumptions in 7.2.1.
\[\frac{dP}{dt} = rP \left( 1 - \frac{P}{K} \right) - \alpha_h P H\] and \[\frac{dH}{dt} = \alpha_h PH - \delta_h H\] Can find two steady states, one with \(H=0\) and \(P = K\). The other is: \[P^* = \delta_h / \alpha_h, \quad H^* = \frac{r}{\alpha_h} \left(1 - \frac{P^*}{K} \right)\]
### params for simulations
a_h <- 0.01 ### herbivore attack rate
r_p <- 0.5 ### plant intrinsic growth rate
d_h <- 0.5 ### carbon mass losses for herbivores
K_p <- 200 ### carrying capacity for plants
pops_0 <- c(P = 25, H = 5) ### initial carbon densities of herbivore and plant
prms <- c(a_h = a_h, r_p = r_p, d_h = d_h, K_p = K_p)
t <- seq(0, 100, by = 1)
foodchain2_out <- ode(pops_0, t, foodchain_2level, prms)
foodchain2_df <- foodchain2_out %>%
as.data.frame() %>%
setNames(c('t', 'P', 'H'))### Equilibrium
P_star2 <- d_h / a_h
H_star2 <- r_p / a_h * (1 - P_star2 / K_p)
### Isoclines
P_iso2 <- function(p) {
h <- r_p / a_h * (1 - p / K_p)
}
H_iso2 <- P_star2
### plot
ggplot(data = foodchain2_df, aes(x = P, y = H)) +
ggtheme_plot() +
stat_function(aes(x = P), fun = P_iso2,
linetype = 'dashed', color = 'darkred') +
geom_vline(xintercept = H_iso2,
linetype = 'dashed', color = 'darkblue') +
geom_segment(aes(xend = lead(P), yend = lead(H)),
arrow = arrow(length = unit(0.07, 'inches')),
size = .25) +
geom_point(x = P_star2, y = H_star2, size = 3) +
geom_hline(yintercept = 0) + geom_vline(xintercept = 0) +
ylim(c(0, max(foodchain2_df$H)))From the equil formulas above, \(P^* = \frac{\delta_h}{\alpha_h}\) = 50, \(H^* = \frac{r}{\alpha_h} \left(1 - \frac{P^*}{K} \right\) = 37.5
To know how each population growth rate changes in response to changes in the abundance each of the other population, we take the partial derivatives.
\[\mathbf J = \begin{pmatrix} \frac{\partial \dot P}{\partial P} & \frac{\partial \dot P}{\partial H} \\ \frac{\partial \dot H}{\partial P} & \frac{\partial \dot H}{\partial H} \end{pmatrix} = \begin{pmatrix} r \left(1 - \frac{2P}{K} \right) - \alpha_h H & -\alpha_h P \\ \alpha_h H & \alpha_h P - \delta_h \end{pmatrix}\]
We can replace the \(P\) and \(H\) in the Jacobian with the equilibria found above \(P^*, H^*\).
\[\begin{pmatrix} r \left(1 - \frac{2\frac{\delta_h}{\alpha_h}}{K} \right) - \alpha_h \frac{r}{\alpha_h} \left(1 - \frac{\frac{\delta_h}{\alpha_h}}{K} \right) & -\alpha_h \frac{\delta_h}{\alpha_h} \\ \alpha_h \frac{r}{\alpha_h} \left(1 - \frac{\frac{\delta_h}{\alpha_h}}{K} \right) & \alpha_h \frac{\delta_h}{\alpha_h} - \delta_h \end{pmatrix} = \begin{pmatrix} \frac{-r \delta_h}{\alpha_h K} & -\delta_h\\ r \left(1 - \frac{\delta_h}{\alpha_h K} \right) & 0 \end{pmatrix}\]
Recall the Routh-Hurwitz criterion:
\[ J_{11} + J_{22} < 0\] \(J_{11}\) is negative, so the plant population will self regulate, while \(J_{22}\) is zero, so the herbivore population does not have a direct density-dependence effect on itself.
\[ J_{11}J_{22} - J_{12}J_{21} > 0\] \(J_{11}J_{22} = 0\), so stability depends on \(J_{12}J_{21} < 0\). The plant declines due to the herbivore (\(J_{12} < 0\)) and the herbivore increases due to the plant (\(J_{21} > 0\) as long as \(P^* < K\)).
\[\begin{pmatrix} \frac{-r \delta_h}{\alpha_h K} & -\delta_h\\ r \left(1 - \frac{\delta_h}{\alpha_h K} \right) & 0 \end{pmatrix}\]
# We can perform eigenanalysis given the parameters above.
Jac <- matrix(c(-r_p * d_h / (a_h * K_p), -d_h, r_p * (1 - d_h / (a_h * K_p)), 0),
byrow = TRUE, nr = 2)
eig <- round(eigen(Jac)[["values"]],3)If we performed eigenanalysis on the above Jacobian matrix, we would find that the eigenvalues are complex (-0.062+0.428i, -0.062-0.428i), meaning that the populations will oscillate or cycle, with period \(\frac{2\pi}{\omega}\). * If the real parts are zero, the system will exhibit neutral stability * If the real parts are negative, the system will converge toward the equilibrium. * If the real parts are positive, the system will approach a stable oscillation away from the equilibrium.
Add a third group, “consumer,” which only eat herbivores, to the ecosystem. Carbon biomass density of consumers is \(C(t)\) gC/m2; respiration & mortality loss at \(\delta_c\) day-1, and linear functional response attack rate \(\alpha_c\) m^2/day/gC. The system of equations:
\[\frac{dP}{dt} = rP \left( 1 - \frac{P}{K} \right) - \alpha_h PH\] \[\frac{dH}{dt} = \alpha_h PH - \delta_h H - \alpha_c H C\] \[\frac{dC}{dt} = \alpha_c HC - \delta_c C\] Can find three steady states:
This will calculate \(\frac{dC}{dt}, \frac{dH}{dt}, \frac{dP}{dt}\).
To have sensible equilibria, the following needs to hold true:
\[K \left(1 - \frac{\alpha_h \delta_c}{\alpha_c r} \right) \geq \frac{\delta_h}{\alpha_h}\]
### params for simulations: what are some sensible numbers?
a_h <- 0.005 ### herbivore attack rate
a_c <- 0.015 ### consumer attack rate
r_p <- 0.5 ### plant intrinsic growth rate
K_p <- 500 ### plant carrying capacity
d_h <- 0.90 ### carbon mass losses for herbivores
d_c <- 0.90 ### carbon mass losses for consumers
K_test <- K_p * (1 - (a_h * d_c) / (a_c * r_p)) >= d_h / a_h
if(!K_test) {
stop('parameters are problematic')
}
### Equilibrium
H_star3 <- d_c / a_c
P_star3 <- K_p * (1 - a_h * H_star3 / r_p)
C_star3 <- (r_p * P_star3 * (1 - P_star3 / K_p) - d_h * H_star3) / d_c
if(any(c(P_star3 < 0, H_star3 < 0, C_star3 < 0))) {
stop('Param values lead to non-sensible equilibria')
}pops_0 <- c(P = 100, H = 100, C = 5)
### initial carbon densities of herbivore, plant, and consumer
prms <- c(a_h = a_h, a_c = a_c, r_p = r_p,
d_h = d_h, K_p = K_p, d_c = d_c)
t <- seq(0, 100, by = 0.1)
foodchain3_out <- ode(pops_0, t, foodchain_3level, prms)
foodchain3_df <- foodchain3_out %>%
as.data.frame() %>%
setNames(c('t', 'P', 'H', 'C'))equil_df <- data.frame(P = P_star3, H = H_star3, C = C_star3) %>%
gather(level, pop) %>%
mutate(level = fct_inorder(level))
pop_v_time_df <- foodchain3_df %>%
gather(level, pop, -t) %>%
mutate(level = fct_inorder(level))
ggplot(pop_v_time_df, aes(x = t, y = pop, color = level)) +
ggtheme_plot() +
geom_line() +
geom_hline(data = equil_df, aes(yintercept = pop, color = level),
linetype = 'dashed', size = .25) +
geom_hline(yintercept = 0) +
scale_color_manual(values = c('darkgreen', 'blue', 'darkred'))Add a fourth group, “top predators,” which only eat consumers, to the ecosystem. Carbon biomass density of top predators is \(T(t)\) gC/m2; respiration & mortality loss at \(\delta_t\) day-1, and linear functional response attack rate \(\alpha_t\) m^2/day/gC. The system of equations:
\[\frac{dP}{dt} = rP \left( 1 - \frac{P}{K} \right) - \alpha_h PH\] \[\frac{dH}{dt} = \alpha_h PH - \delta_h H - \alpha_c H C\] \[\frac{dC}{dt} = \alpha_c HC - \delta_c C - \alpha_c CT\] \[\frac{dT}{dt} = \alpha_t CT - \delta_t T\] Can find four steady states:
This will calculate \(\frac{dC}{dt}, \frac{dH}{dt}, \frac{dP}{dt}\).
foodchain_4level <- function(t, y, params) {
P <- y[1]
H <- y[2]
C <- y[3]
TT <- y[4]
with(as.list(params), {
dP_dt <- r_p * P * (1 - P / K_p) - a_h * H * P
dH_dt <- a_h * P * H - d_h * H - a_c * H * C
dC_dt <- a_c * H * C - d_c * C - a_t * C * TT
dT_dt <- a_t * C * TT - d_t * TT
return(list(c(dP_dt, dH_dt, dC_dt, dT_dt)))
})
}To have sensible equilibria, the following needs to hold true: ??? update \[\phi \geq \delta_p \left(\frac{\delta_h}{\alpha_h}\right) + \delta_h \left(\frac{\delta_c}{\alpha_c} \right)\]
### params for simulations: what are some sensible numbers?
a_h <- 0.015 ### herbivore attack rate
a_c <- 0.025 ### consumer attack rate
a_t <- 0.035 ### top predator attack rate
r_p <- 0.25 ### intrinsic growth rate for plants
d_h <- 0.90 ### carbon mass losses for herbivores
K_p <- 500 ### carrying capacity for plants
d_c <- 0.90 ### carbon mass losses for consumers
d_t <- 0.90 ### carbon mass losses for top predators
### Equilibrium
C_star4 <- d_t / a_t
P_star4 <- (d_h + a_c * C_star4) / a_h
H_star4 <- r_p / a_h * (1 - P_star4 / K_p)
T_star4 <- (a_c * H_star4 - d_c) / (a_t)
# if(any(c(P_star4 < 0, H_star4 < 0, C_star4 < 0, T_star4 < 0))) {
# stop('Param values lead to non-sensible equilibria')
# }pops_0 <- c(P = 100, H = 100, C = 20, TT = 10)
### initial carbon densities of herbivore, plant, and consumer
prms4 <- c(a_h = a_h, a_c = a_c, a_t = a_t, r_p = r_p,
d_h = d_h, K_p = K_p, d_c = d_c, d_t = d_t)
t <- seq(0, 100, by = 0.1)
foodchain4_out <- ode(pops_0, t, foodchain_4level, prms4)
foodchain4_df <- foodchain4_out %>%
as.data.frame() %>%
setNames(c('t', 'P', 'H', 'C', 'T'))equil_df <- data.frame(P = P_star4, H = H_star4, C = C_star4, T = T_star4) %>%
gather(level, pop) %>%
mutate(level = fct_inorder(level))
pop_v_time_df <- foodchain4_df %>%
gather(level, pop, -t) %>%
mutate(level = fct_inorder(level))
ggplot(pop_v_time_df, aes(x = t, y = pop, color = level)) +
ggtheme_plot() +
geom_line() +
geom_hline(data = equil_df, aes(yintercept = pop, color = level),
linetype = 'dashed', size = .25) +
geom_hline(yintercept = 0) +
scale_color_manual(values = c('darkgreen', 'blue', 'red', 'darkred'))Plant carbon biomass as constant primary production (same as 7.2.1)
\[\frac{dP}{dt} = \phi - \delta_p P\]
Add a second group, “herbivores,” to the ecosystem, but now the herbivores have a type II functional response.
\[\frac{dP}{dt} = \phi - \delta_p P - \frac{\alpha_h H P}{P + P_0}\] and \[\frac{dH}{dt} = \frac{\alpha_h PH}{P + P_0} - \delta_h H\] Can find two steady states, one with \(H=0\) and \(P = \phi/\delta_p\). The other is:
\[P^* = \frac{\delta_h P_0}{\alpha_h - \delta_h}, \quad H^* = \frac{\phi - \delta_p P^*}{\delta_h}\]
These should be positive, as long as:
\[\frac{\phi}{\delta_p} \geq \frac{\delta_h P_0}{\alpha_h - \delta_h}\]
### params for simulations
a_h <- 0.75 ### herbivore attack rate
phi <- 150 ### plant constant primary production
d_h <- 0.32 ### carbon mass losses for herbivores
d_p <- 0.25 ### carrying capacity for plants
P_0 <- 200
pops_0 <- c(P = 25, H = 5) ### initial carbon densities of herbivore and plant
prms <- c(a_h = a_h, phi = phi, d_h = d_h, d_p = d_p, P_0 = P_0)
t <- seq(0, 100, by = 1)
foodchain2_out <- ode(pops_0, t, foodchain_2level, prms)
foodchain2_df <- foodchain2_out %>%
as.data.frame() %>%
setNames(c('t', 'P', 'H'))### Equilibrium
P_star2 <- d_h * P_0 / (a_h - d_h)
H_star2 <- (phi - d_p * P_star2) / d_h
### Isoclines
P_iso2 <- function(p) {
h <- (phi - d_p * p) / (a_h * p) * (p + P_0)
}
H_iso2 <- d_h / a_h * P_0 / (1 - d_h / a_h)
### plot
ggplot(data = foodchain2_df, aes(x = P, y = H)) +
ggtheme_plot() +
stat_function(aes(x = P), fun = P_iso2,
linetype = 'dashed', color = 'darkred') +
geom_vline(xintercept = H_iso2,
linetype = 'dashed', color = 'darkblue') +
geom_segment(aes(xend = lead(P), yend = lead(H)),
arrow = arrow(length = unit(0.07, 'inches')),
size = .25) +
geom_point(x = P_star2, y = H_star2, size = 3) +
geom_hline(yintercept = 0) + geom_vline(xintercept = 0) +
ylim(c(0, max(foodchain2_df$H)))The “consumer” functional group now has a type II functional response as well.
As usual, can find three steady states; coexistence steady state, with
This will calculate \(\frac{dC}{dt}, \frac{dH}{dt}, \frac{dP}{dt}\).
foodchain_3level <- function(t, y, params) {
P <- y[1]
H <- y[2]
C <- y[3]
with(as.list(params), {
dP_dt <- phi - d_p * P - (a_h * H * P) / (P + P_0)
dH_dt <- (a_h * H * P) / (P + P_0) - d_h * H - (a_c * C * H) / (H + H_0)
dC_dt <- (a_c * C * H) / (H + H_0) - d_c * C
return(list(c(dP_dt, dH_dt, dC_dt)))
})
}To have sensible equilibria, the following needs to hold true:
\[K \left(1 - \frac{\alpha_h \delta_c}{\alpha_c r} \right) \geq \frac{\delta_h}{\alpha_h}\]
### params for simulations: what are some sensible numbers?
a_h <- 0.75 ### herbivore attack rate
a_c <- 0.4 ### consumer attack rate
phi <- 150 ### plant constant primary production
d_h <- 0.32 ### carbon mass losses for herbivores
d_p <- 0.25 ### carrying capacity for plants
P_0 <- 200
H_0 <- 150
d_c <- 0.20 ### carbon mass losses for consumers
### Equilibrium
H_star3 <- (d_c * H_0) / (a_c - d_c)
### params for quadratic formula:
a <- d_p
b <- - (phi - d_p * P_0 - a_h * H_star3)
c <- -phi * P_0
P_star3 <- (-b + sqrt(b^2 - 4 * a * c)) / (2 * a)
C_star3 <- (phi - d_p * P_star3 - d_h * H_star3) / d_c
if(any(c(P_star3 < 0, H_star3 < 0, C_star3 < 0))) {
stop('Param values lead to non-sensible equilibria')
}pops_0 <- c(P = 100, H = 100, C = 5)
### initial carbon densities of herbivore, plant, and consumer
prms <- c(a_h = a_h, a_c = a_c, phi = phi,
d_h = d_h, d_p = d_p, d_c = d_c, P_0 = P_0, H_0 = H_0)
t <- seq(0, 250, by = 0.1)
foodchain3_out <- ode(pops_0, t, foodchain_3level, prms)
foodchain3_df <- foodchain3_out %>%
as.data.frame() %>%
setNames(c('t', 'P', 'H', 'C'))equil_df <- data.frame(P = P_star3, H = H_star3, C = C_star3) %>%
gather(level, pop) %>%
mutate(level = fct_inorder(level))
pop_v_time_df <- foodchain3_df %>%
gather(level, pop, -t) %>%
mutate(level = fct_inorder(level))
ggplot(pop_v_time_df, aes(x = t, y = pop, color = level)) +
ggtheme_plot() +
geom_line() +
geom_hline(data = equil_df, aes(yintercept = pop, color = level),
linetype = 'dashed', size = .25) +
geom_hline(yintercept = 0) +
scale_color_manual(values = c('darkgreen', 'blue', 'darkred'))